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A simple fragmentation model is introduced and analysed. We show that, under very general 
conditions, an effective power law for the mass distribution arises with realistic exponent. This ex- 
ponent has a universal limit, but in practice the effective exponent depends on the detailed breaking 
mechanism and the initial conditions. This dependence is in good agreement with experimental 
results of fragmentation. 
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One of the best known physical processes in Nature is fragmentation. From our daily experience we know that a 
^\ \ material bulk under stress or shock will break into smaller pieces. Experiments show generally that the number of 
t-H fragments with a linear size larger than r behaves like 



=5 



O 

(N 



O 



N{r)~r- U . (1) 



The exponent D, usually called fractal dimension, is found to lie in the range D m 2 -f- 3 for fragmentation of 3- 
dimensional objects Q. A theoretical understanding of the statistical origin of eq. (|l|) is currently pursuited by many 
authors. The simplest models || predict a log-normal distribution, incompatible with observation. More refined 
I | models with various assumptions about breaking mechanism yield a power law behavior for N(r). However 

^ . it appears that theories producing a single universal exponent may not account for the experimental range of D. 
ON ' Recently quite large amount of numerical simulations Q with rather realistic physical parameters such as stress, 
shear, neighborhood were able to show qualitatively correct power laws. 

In this Letter we develop an analytical model of fragmentation, without using specific breaking mechanism. We 
shall see that under very general, simple conditions an effective power law will result. A general theory is called for, 
since it is difficult to observe what is really happening during fragmentation. We imagine that fragmentation happens 
in a hierarchical order — i.e. a large piece breaks into n smaller ones and each of the n pieces may then break further. 
. For simplicity we assume n = 2 for all d and levels. Our numerical simulations show that other finite n does not affect 
d ' our conclusions. 

At level k of the hierarchy, we consider one object of volume V and energy E, the only variables retained in our 
model. E is the total energy, including kinetic energy, elastic energy, etc. We stress, however, that at our level of 
description this specification is not necessary. The only property of the energy we retain is that it is conserved. That 
t£ is not exactly true in real fragmentation but it is a good approximation for large masses because usually dissipation is 
an effect proportional to the area of the fracture produced. We assume further that it is the energy density E/V that 
decides whether an object breaks: if E/V exceeds a threshold, set to 1 for all k, the object breaks further, otherwise 
not. If E/V > 1, the object breaks into two pieces of energy e and E — e and volume v and V — v respectively. The 
two resulting fragments will or will not break in their turn according to their energy-volume ratios, i.e. if 

x x = - > 1 ; x 2 = §^1 > 1. (2) 
v V — v 

The above process is repeated for an arbitrary number k of levels. Note that k is not necessarily proportional to time. 
The above variables are all for level k and we have suppressed the sub-index k for clarity. At level k — 0, E and Vo 
are given by the initial energy and volume of the system (Eq/Vq > 1). 

Let q(e, v\E)dedv be the probability that the energy and volume of an element are between e and e + de and v and 
v + dv, given that it results from the fragmentation of an object of unit volume V = 1 and energy E > 1 (for our 
purposes we may consider the volume V to be unity at any level, since only the ratio E/V matters). This distribution 
accounts for all the informations of a detailed breaking mechanism. We shall assume that q(e, v\E) — -^ip(v) ||. 
This implies an uniform distribution in energy, but arbitrary distribution in volume (for symmetry we require (p(v) = 
ip(l — v)). This certainly is not the most general case, but it still includes a large class of models. 

The fragments with energy density x > 1 are called "unstable" , those with x < 1 are called "stable" . The nature of 
the cascade process is best illustrated by studying the distribution pk (x) of energy density x of the unstable elements 
at the level k. This can be computed from the initial distribution po(x) once a recurrence relation between Pk(x) and 
Pk+i(x) is found. In order to derive this relation, let us consider a unit volume object with energy E > 1 breaking in 
two smaller elements. The jacobian of the transformation (e, v) — ► {x\ 1 xi) is readily found using (2), and we find the 
joint distribution of x\ and x 2 from the above q(e,v\E), 
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Denote x+ = maxfii, x 2 ) and X— = min(xi, X2), from (2) we have X- < E < x + . Integrating out one of the two 
variables we obtain the distributions of x±, respectively: 

, . , 20(x - E) f E/x , . , 
p+[x\E) = / vtp{v)dv 

P-{x\E) = - (4) 

where 6(x) — for x < and 9(x) = 1 otherwise. The + element always breaks further since x+ > E > 1, while the 
— element does not break if X- < 1, this occurs with the probability 



C{E) = J p_(x\E)dx = ± 



(5) 



Let us define the distribution of the energy density p(x\E) of one unstable (x > 1) element between the two 
fragments. The pair x+ falls into two cases: i) If x_ < 1, which occur with probability C(E) given above, the 
only possibility is that x — x + . In this case the distribution p(x\E) is just that of x + , i.e. p + (x\E). ii) If x- > 1, 
which occur with probability 1 — C{E) , x can be either x+ or In this case, occurring with the probability 1 — C [E) , 
p(x\E) is the average of the distributions p±(x\E) with weight 1/2 each. Taking into account the two cases we find: 

E + 1 f E / x 1 
p(x\E) = 9{x - E)—^— / vcp(v)dv + 6{E - x)- 



One can verify that p(x\E) is normalized in [1, 00). Using p(x\E) we finally construct the iteration relation: 



p k+ i{x) = J dEp(x\E)p k (E). (6) 
As k — > 00 the distribution pk(x) is expected to converge to a limit which describes the asymptotic behavior: 

/>oo 

Poo(x) = / p(x\E) PoD (E)dE. (7) 



The solution of this equation will depend on the specific choice of <p(v). For the special case <p(v) = 1, i.e. the 
uniform distribution which is interesting per se, eq. (0) can be solved to give: Poo(x) — ^ (l — ^) exp (— ^-). This 
distribution is not normalizable because Poo(x) ~ A/x for large x. For arbitrary ^(f), detailed study of eq. (^) for 
x — > 00 reveals that this asymptotic behavior holds as well. Indeed setting p x (x) ~ Ax~ 7 in eq. (Q) and carrying out 
an asymptotic analysis, one finds that the exponent 7 must satisfy the equation 

7= f v^if^dv. (8) 
Jo 

This equation has always a solution for 7=1. For convex distributions <p(v), i.e. they weight more the center of 
the interval [0, 1], there is a second solution 7 > 1. However the solution with 7=1 clearly dominates the large x 
behavior. For concave <p(v) the second solution becomes 7 < 1. We expect the physical relevant distributions are 
convex since objects tend to break in the middle easier than at the edges. 

Denote by Ck the probability that one of the fragments (the other is unstable by definition) becomes stable at level 

k; 

/CO 
C(E) Pk (E)dE. (9) 

The important consequence of 7 = 1 is that Ck — > as k — > 00, regardless of the distribution <p{v). To see this, we 
assume for simplicity that Pk(x) — A^x tk with 7^ — > 1 + as k — > 00. Then the integral in (Q) is finite for all 7fe > 1. 
On the other hand, as 7^. — » 1 + , the constant which enforces normalization vanishes. Therefore Ck ~ vanishes 
as well. 
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Unfortunately we cannot obtain C k in closed form even for the uniform distribution ip — 1. The numerical procedure 
is straightforward: start with a distribution at level k = with an initial condition po(E) = 8{E — Eq), Eq > 1, iterate 
eq. (6) to the desired level k, then find Ck using (9). The numerical results for a particular ip and various initial 
energies are plotted in Fig. [j]. The function Ck is not universal. It depends on Eq as well as on the function ip. In 
general, for a given ip, Ck(Eo) decreases with Eo for fixed k, and vanishes slowly when k — > oo (see fig. ^). 

From the above preparations, we can establish the desired scaling laws. In eq. (|l|) only the stable fragments left in 
the cascade count. N(r) is trivially related to the distribution W(V) of volumes V = r d : N(r) — J v W(V)dV where 
W(V)dV is the number of stable fragments with volume between V and V + dV . We thus have the relation from (1): 

W(V) oc V- a -\ D = da. (10) 

W(V) receives contributions Wk(V) from all fragmentation levels and it can be written: 

oo oo 

W(V) =J2 w k(V) ^J2 C k N kW k (V), (11) 
fc=l fc=l 

Here Wk(V) is the volume distribution of all the fragments (stable and unstable) produced by the fc th step of frag- 
mentation, Nk is the number of unstable fragments at level k. Nk unstable objects produce 2Nk fragments, of which 
(2 — Ck)Nk 1 which is Nk+i, are unstable again. Thus CkNk stable fragments are produced at level k. Eq. ( |Tl| ) 
assumes that stable fragments are produced at level k with a probability Ck which is independent of their volume 
V. This is clearly an approximation because in reality Ck — Ck(V) Neglecting this dependence allows us to keep 
the discussion at an elementary level. Moreover for distributions tp(v) which are sharply peaked around v = 1/2 this 
approximation is reliable. Indeed for tp(v) = 5(v — 1/2) the dependence on the volume clearly drops (all the fragments 
at level k have the same volume). For general distributions this approximation yields an upper bound to the true 
exponent a ||. We shall see below that even for broad distributions, such as the uniform one, we recover qualitatively 
correct results. As we shall see the fractal dimension is determined by the exponential behavior of Nk with k, which 
depends on Ck, and the explicit factor Ck in eq. (^) does dot play any role. This observation supports the present 
approximation. Moreover the same qualitative behavior discussed below results from a more detailed calculation, to 
be presented elsewhere 0] . 

For the uniform case (Mv) = 1 it is easy to find [|| that w k (V) = (\nV) k ^ 1 /(k - 1)!. Assuming C k N k oc (2 - C*) k , 
one can easily sum eq. ( jll| ) with the result W(V) oc V~ a ~ x with a = 1 — C* . Let us generalize this analysis for an 
arbitrary distribution ip(v). Let us evaluate the m th moment of V using (11) and (J10|) . Note that some moments 



can be divergent and there is a smallest to* for this to happen. Multiply both sides of (|llj) by V m and integrate 
over V. On the RHS one finds (V m )k where the average is done using the distribution Wk(V). Since V, within our 

approximation, is the product of k independent variables, each distributed by ip, we have (V m )k = [J v m ip{v)dv\ k . 
Multiplying (V m )k by CkNk and summing over k, we see (e.g. by the ratio method) that the sum first diverges when 
(2 — C*) J v m " tp(v)dv — 1, where C* is defined by N k +i = (2 - C*)N k - This divergence must be matched to the one 
occurring on the LHS of the equation, which is proportional to J dVV™ 1 ^ " 1 . The latter occurs for m* = a, which 
therefore gives the relation 



v a p(v)dv = j—^ 7 (12) 



which implicitly yields the exponent a within the present approximation. For <p = 1, eq. (|T^) reduces to the result 
we found previously. Most importantly, for C* — 0, and only for this value, one finds a = 1, regardless of ip\ 

In order to derive the value of C* , which determines the scaling exponent, we observe that in real life the observation 
of W(V) is limited to a finite window V m i n < V < Vo, which may cover several decades. C k is a slowly varying function 
of k whereas the typical V decreases exponentially with k. Thus, over an exponentially large range of V, Ck can 
be regarded as constant and an effective power law can be established. More precisely, W / (V r ) in this window, is 
dominated by contributions around a certain k* (see fig. |l|) thus giving C* = Ck'- Further fragmentation for k > k* , 
seldom adds stable fragments larger than V m i n . Pursuing this arguments further, one can find that, to leading order, 
k* ~ clog(Vo/V^ n i n ) where c is a coefficient of order unity (c = 1 + a for the uniform case within our approximation). 
k* appears to depend weakly on the initial energy Eq, as shown in Fig.l. Therefore, strictly speaking, there is no true 
power law, except at the ideal limit k* — ► oo where one can take C* = so that D = d. 

How does the exponent a vary in practical situations? First we notice that often, in experiments, Kni n is set by 
a dissipation scale, below which an object cannot break anymore. In other words at this scale the energy lost in 
the rupture of an object becomes non-negligible with respect to its energy. Thus we conclude that a depends on 
Vq through the dependence of C* = Ck* on k* ~ log(Vb/Kni n ). This suggests that for Vo/V m i n — > oo it is possible 
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to recover the ideal limit D = d. Moreover C* depends also on E$. A larger Eq results in a smaller C* , since Cfc 
decreases with Eq (see fig. 0). Thus the effective a(Eo) increases with Eq, until it reaches the universal value of 1. 

Figure || shows that, in spite of the approximation used to derive eq. |l2L this scenario is confirmed by numerical 
simulations. It also has many features that have been observed in real experiments. In fragmentation of 2 dimensional 
glass plates Q , the fractal dimension was nicely extrapolated to D = 2 for infinite input energy The same agreement 
can be found comparing our results with 3 dimensional data jj]. We note in particular that D for stony meteorites 
and asteroids, which are the result of an extremely long fragmentation process, where very high energies are involved, 
are in excellent agreement with D = 3. For these systems we expect Vo/Vmm and Eq/Vo to be extremely large. On 
the other hand, in projectile fragmentation the input energy is the kinetic energy of the projectile. The volume to 
be considered is the total volume of the system, which is essentially that of the target object. Therefore, even for 
extremely high kinetic energies, the input energy density E /Vq which would enter our calculation can be relatively 
small. Indeed, experiments of projectile fragmentation yield exponents D rs 2.5, below the universal value D = 3. 
Our analysis yields d as an upper bound for D (C* > implies a < 1). It is reassuring that, apart from the 
data concerning ash and pumice, materials which probably need a separate treatment, all the data in rcf. jjj have 
D < d = 3. Comparing the exponents obtained for different materials, in terms of our model, translates into computing 
the exponent for different distributions (p( v )- As shown in Fig. ^|, for convex distributions the convergence of a(E) 
to 1 is faster than that obtained for concave distributions. This might be related to the fact that, from eq. (|12j), 
a ~ 1 — C* /\A(vlnv) \ + 0(C* 2 ), where the coefficient of C* is larger for broader distributions. 

Problems arise when considering the d — 1 experiments described e.g. in ref. [Io| . There it was found that long 
thin glass rods fragmentation produces a size distribution with an exponent D ss 1.5. This is clearly inconsistent with 
our analysis since it would need a C* < 0. This failure, we believe, lies in the assumption that the breaking of a large 
object is determined by its global energy density. Without loss of generality we know that energy correlation inside 
a volume propagates via nearest neighbor interaction, this leads to a Laplace equation. For correlation induced by a 
Laplace equation we know that for d > 2 the correlation is very strong and d = 2 is the marginal case. This strong 
correlation allows to describe the fragmentation of one object as an event which produces two objects and which 
depends on a single variable, its energy density. For d — 1 this is not true. The energy is very loosely correlated along 
a line. We suspect that, because of the weak correlation, simultaneous breaking will happen in many uncorrclatcd 
regions of a large d = 1 object, making our scenario invalid. For smaller and smaller rod length, the energy correlation 
becomes stronger and stronger. We therefore expect that below a certain length threshold, our scenario can be applied. 
Remarkably, in experiments of fragmentation of long glass rods ||l0[| , a crossover occurs and for intermediate sizes the 
mass distribution is described by an exponent D « 0.6 < 1. 

In this work we have analysed a simple fragmentation model. We show that under very general conditions an effective 
power law arises. The exponent is not universal but depends on the detailed mechanism and the initial conditions. 
There is an ideal universal limit, independent of any of our choices, which can be approached for higher input energies. 
This qualitative behavior, derived within a simple approximation for the mass distribution, is completely consistent 
with numerical simulations of our model. Furthermore our results agree reasonably well with existing experiments in 
2-d and 3-d. We also argue that our model is not directly applicable in d = 1. A more precise quantitative theory, 
assuming a general distribution q(s,v), going beyond the approximations in eq. (11) and including explicitly the 
effects of dissipation, is under current investigation 0j . 
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FIG. 1. Cfe as a function of k for Vb = 1 and Eo = 2 (O), Eo — 8 (□) and So = 32 (o). The data refers to the tent distribution 
<p(v) — 2v for v < 1/2 and tp(l — v) = f(v). The lines refers to the contribution J y Wk(V)dV of the level k to the statistics 

of W(V) for Vmin < V < 1 with Knin = 3 ■ 10~ 5 . These lines refer, from left to right, to the same values of Eo — 2, 8 and 32 
used for Ck- 

FIG. 2. Fit of the mass size distribution for <p(v) = 1 and Eo/Vo = 4 (□) and for the tent distribution and Eq/Vo = 8 (O) 
and for ip(v) = S(v — 1/2) and Eo = 16 (o). Inset: exponent a, obtained from the fit, versus Eo/Vo for the uniform (□), the 
tent (O) and the delta (o) distributions. 
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